! (c) British Crown Copyright, the Met Office.
! All rights reserved.
! 
! Redistribution and use in source and binary forms, with or without modification, are permitted 
! provided that the following conditions are met:
! 
!     * Redistributions of source code must retain the above copyright notice, this list 
!       of conditions and the following disclaimer.
!     * Redistributions in binary form must reproduce the above copyright notice, this list
!       of conditions and the following disclaimer in the documentation and/or other materials 
!       provided with the distribution.
!     * Neither the name of the Met Office nor the names of its contributors may be used 
!       to endorse or promote products derived from this software without specific prior written 
!       permission.
! 
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


MODULE MOD_COSP_RTTOV
  USE rttov_const, ONLY : errorstatus_success, errorstatus_fatal

  ! rttov_types contains definitions of all RTTOV data types
  USE rttov_types, ONLY :     &
  rttov_options,       &
  rttov_coefs,         &
  profile_type,        &
  transmission_type,   &
  radiance_type,       &
  rttov_chanprof,      &
  rttov_emissivity

  USE rttov_unix_env, ONLY : rttov_exit

  USE PARKIND1, ONLY : jpim, jprb, jplm

  IMPLICIT NONE

  INCLUDE "rttov_direct.interface"
  INCLUDE "rttov_alloc_prof.interface"
  INCLUDE "rttov_alloc_rad.interface"
  INCLUDE "rttov_alloc_transmission.interface"
  INCLUDE "rttov_dealloc_coefs.interface"
  INCLUDE "rttov_user_options_checkinput.interface" 
  INCLUDE "rttov_read_coefs.interface" 

  INTEGER, PARAMETER :: MaxLim  =  10000

CONTAINS
  SUBROUTINE RTTOV_MULTPROF( &
        nch_in,     & ! number of channels 
        ichan_in,   & ! channel indices
        surfem_in,  & ! surface emissivity values 
        prf_num_in, & ! number of profiles to simulate
        nlevels_in, & ! number of pressure levels
        plat_in,    & ! platform number
        sat_in,     & ! satellite number
        sens_in,    & ! instrument number
        zenang_in,  & ! zenith angle 
        p_in,       & ! Pressure [hPa]
        t_in,       & ! Temperature [ K ]
        q_in,       & ! Specific humidity [ ppmv ]
        o3_in,      & ! Ozone VMR [ ppmv ]
        co2_in,     & ! CO2 VMR [ ppmv ] *This is a single value* 
        ch4_in,     & ! CH4 VMR [ ppmv ] *This is a single value*
        n2o_in,     & ! N2O VMR [ ppmv ] *This is a single value*
        co_in,      & ! CO VMR [ ppmv ]  *This is a single value*
        h_surf,     & ! Surface height [ m ]
        u_surf,     & ! U wind at 10 m  [ m/s ]
        v_surf,     & ! V wind at 10 m [ m/s ]
        t_skin,     & ! Skin Temperatre [ K ] 
        p_surf,     & ! Surface pressure 
        t_surf,     & ! 1.5 m Temperature [ K ]
        q_surf,     & ! 1.5 m specific humidity [ ppmv ]
        lsmask,     & ! land sea mask
        latitude,   & ! Latitude [ deg north ]
        tbs         & ! Brightness temperature [ K ] (OUTPUT)
        )
    
    !------ Input arguments. No RTTOV kinds should be used here -----------------
    INTEGER, INTENT(in)  :: nch_in             ! number of channels to be computed
    INTEGER, INTENT(in)  :: ichan_in(nch_in)   ! Indices of selected channels 
    REAL, INTENT(in)     :: surfem_in(nch_in)  ! Surface emissivities for the channels
    INTEGER, INTENT(in)  :: prf_num_in
    INTEGER, INTENT(in)  :: nlevels_in
    INTEGER, INTENT(in)  :: plat_in  ! Satellite platform
    INTEGER, INTENT(in)  :: sat_in   ! Satellite number
    INTEGER, INTENT(in)  :: sens_in  ! Satellite sensor
    REAL, INTENT(in)     :: zenang_in          ! Satellite zenith angle

! DAR: Probably sub-optimal to allocate automatic arrays for intent(in) variables?
!      Prefer to use assumed size arrays from calling subroutine but will leave for now
    REAL, INTENT(in)     :: p_in(prf_num_in, nlevels_in)  ! Pressure profiles  
    REAL, INTENT(in)     :: t_in(prf_num_in, nlevels_in)  ! Temperature profiles
    REAL, INTENT(in)     :: q_in(prf_num_in, nlevels_in)  ! Humidity profiles
    REAL, INTENT(in)     :: o3_in(prf_num_in, nlevels_in) ! Ozone profiles
!    REAL, INTENT(in)     :: p_in(1:,:)  ! Pressure profiles  
!    REAL, INTENT(in)     :: t_in(1:,:)  ! Temperature profiles
!    REAL, INTENT(in)     :: q_in(1:,:)  ! Humidity profiles
!    REAL, INTENT(in)     :: o3_in(1:,:) ! Ozone profiles

    ! The following trace gases contain constant values
    REAL, INTENT(in) ::  co2_in ! Carbon dioxide 
    REAL, INTENT(in) ::  ch4_in ! methane 
    REAL, INTENT(in) ::  n2o_in ! n2o 
    REAL, INTENT(in) ::  co_in  ! carbon monoxide
    REAL, INTENT(in) ::  h_surf(prf_num_in)         ! Surface height
    REAL, INTENT(in) ::  u_surf(prf_num_in)         ! U component of surface wind
    REAL, INTENT(in) ::  v_surf(prf_num_in)         ! V component of surface wind
    REAL, INTENT(in) ::  t_skin(prf_num_in)         ! Surface skin temperature
    REAL, INTENT(in) ::  p_surf(prf_num_in)                  ! Surface pressure
    REAL, INTENT(in) ::  t_surf(prf_num_in)         ! 1.5 m Temperature
    REAL, INTENT(in) ::  q_surf(prf_num_in)         ! 1.5 m Specific humidity
    REAL, INTENT(in) :: lsmask(prf_num_in)         ! land-sea mask
    REAL, INTENT(in) :: latitude(prf_num_in)       ! Latitude

    REAL, INTENT(inout) :: tbs(prf_num_in, nch_in)  ! Tbs (in the right format)
    
    !------ Local variables. Use only RTTOV kinds or derived types.
    !       logical variables are declared with the same kind
    !       as integers, as they are affected inthe same way by flags like -qintsize=8

    TYPE(rttov_options)                  :: opts           ! Options structure
    TYPE(rttov_coefs),       ALLOCATABLE :: coefs(:)       ! Coefficients structure
    TYPE(rttov_chanprof),    ALLOCATABLE :: chanprof(:)    ! Input channel/profile list
    TYPE(profile_type),      ALLOCATABLE :: profiles(:)    ! Input profiles
    LOGICAL(KIND=jplm),      ALLOCATABLE :: calcemis(:)    ! Flag to indicate calculation of emissivity within RTTOV
    TYPE(rttov_emissivity),  ALLOCATABLE :: emissivity(:)  ! Input/output surface emissivity
    TYPE(transmission_type)              :: transmission   ! Output transmittances
    TYPE(radiance_type)                  :: radiance       ! Output radiances
   
    INTEGER(Kind=jpim), ALLOCATABLE :: instrument(:,:) ! instrument id (3 x n_instruments)
    INTEGER(Kind=jpim), ALLOCATABLE :: nchan(:) ! number of channels per instrument 
    INTEGER(Kind=jpim), ALLOCATABLE :: ichan(:,:)   ! channel list per instrument

    INTEGER(Kind=jpim) :: nprofiles, asw 
    INTEGER(Kind=jpim), ALLOCATABLE :: rttov_errorstatus(:)  ! rttov error return code

    INTEGER(Kind=jpim) :: mxchn 
    INTEGER(Kind=jpim) :: Err_Unit        ! Logical error unit (<0 for default)
    INTEGER(Kind=jpim) :: verbosity_level ! (<0 for default)
    INTEGER(Kind=jpim) :: nrttovid     ! maximum number of instruments
    INTEGER(Kind=jpim) :: no_id        ! instrument loop index
    INTEGER(Kind=jpim) :: i, j, jch
    INTEGER(Kind=jpim) :: Nprof  ! Number of calls to RTTOV
    INTEGER(Kind=jpim) :: nch ! intermediate variable
    INTEGER(Kind=jpim) :: errorstatus, io_status
    INTEGER(Kind=jpim) :: ioout, interp
    INTEGER(Kind=jpim) :: ich, ich_temp, nchanprof, nchannels, prof, chan, oldprof
    INTEGER(Kind=jpim) :: alloc_status(60)
    
!    REAL(Kind=jprb),    ALLOCATABLE :: surfem(:,:),emis(:)   ! surface input emissivity per channel , instrument
!    REAL(Kind=jprb),    ALLOCATABLE :: emissivity (:)
    REAL(Kind=jprb),    ALLOCATABLE :: input_emissivity (:)
   
    CHARACTER (len=14)  :: NameOfRoutine = 'RTTOV_MULTPROF'
    
    LOGICAL :: refrac, solrad, laerosl, lclouds, lsun, all_channels
    
    ! Local variables for input arguments that need type casting to avoid type-mismatch with 
    ! RTTOV kinds. This happens with some compiler flags (-qintsize=8).
    INTEGER(Kind=jpim)  :: prof_num
    INTEGER(Kind=jpim)  :: nlevels

  ! --------------------------------------------------------------------------
  ! 0. Initialise COSP-specific things
  ! --------------------------------------------------------------------------

    ! Type-casting of input arguments that need to be passed to RTTOV
    prof_num = prf_num_in
    nlevels  = nlevels_in
    nprof = prof_num
                
    ! Unit numbers for input/output
    IOOUT = 2
    
    ! Currently we plan to calculate only 1 instrument per call
    nrttovid  =  1    
    mxchn  =  nch_in
    
    errorstatus     = 0
    alloc_status(:) = 0

    ALLOCATE(coefs(nrttovid), stat = alloc_status(1))
    
    ALLOCATE(instrument(3, nrttovid), stat = alloc_status(4))
    
    !maximum number of channels allowed for one instrument is mxchn
!    ALLOCATE(surfem(nch_in, nrttovid), stat = alloc_status(11))
    ALLOCATE(ichan(nch_in, nrttovid), stat = alloc_status(12))

    !DAR FIXME: ichan_in needs an extra dimension so we can determine which
    !           channels belong to which instrument

    DO no_id = 1, NRTTOVID
      ichan(:, no_id)   = ichan_in
!      surfem(:, no_id)  = surfem_in
    ENDDO
   
! --------------------------------------------------------------------------
! 1. Populate RTTOV options structure
! --------------------------------------------------------------------------
    
    opts%config%do_checkinput = .FALSE.
!   opts%rt_ir%addsolar = .TRUE.          ! Include solar radiation
    opts%rt_ir%addsolar = .FALSE.          ! Not supported at the moment

    opts%interpolation%addinterp  = .TRUE.  ! Allow interpolation of input profile
    opts%rt_all%addrefrac         = .TRUE.  ! Include refraction in path calc

    opts%rt_ir%addclouds          = .FALSE. ! Don't include cloud effects
    opts%rt_ir%addaerosl          = .FALSE. ! Don't include aerosol effects
    
!   opts%rt_ir%ozone_data         = .TRUE. ! We have an ozone profile
    opts%rt_ir%ozone_data         = .FALSE.! We have an ozone profile
!   opts%rt_ir%co2_data           = .TRUE. ! We do not have profiles
!   opts%rt_ir%n2o_data           = .TRUE. ! for any other constituents
!   opts%rt_ir%ch4_data           = .TRUE. !
!   opts%rt_ir%co_data            = .TRUE. !
    
    opts%rt_all%use_q2m           = .FALSE.
    opts%rt_mw%clw_data           = .FALSE. !
    
! END OF 1.

    ! DAR: maintain option to use multiple instruments but this will require work on upstream interface
    ! Forward model only

    DO NO_ID = 1, NRTTOVID

  ! --------------------------------------------------------------------------
  ! 2. Read coefficients
  ! --------------------------------------------------------------------------
      
      instrument(1, no_id) = plat_in
      instrument(2, no_id) = sat_in
      instrument(3, no_id) = sens_in

      CALL rttov_read_coefs(errorstatus, coefs(no_id), opts, instrument = instrument(:, no_id))
      CALL rttov_error('fatal error reading coefficients')
            
      ! Ensure the options and coefficients are consistent
      CALL rttov_user_options_checkinput(errorstatus, opts, coefs(no_id))
      CALL rttov_error('fatal error in rttov options')

      IF(ANY(coefs(no_id)%coef%ff_val_chn(:) /= 1 )) THEN
        WRITE(*,*) 'some requested channels have bad validity parameter'
        
        DO ich = 1, SIZE(coefs(no_id)%coef%ff_val_chn(:))
          WRITE(*,*) ich, coefs(no_id)%coef%ff_val_chn(ich)
        END DO
      ENDIF
    ENDDO

! End of 2.

! --------------------------------------------------------------------------
! 4. Allocate profiles, radiance and transmittance structures
! --------------------------------------------------------------------------
        
    asw = 1 ! Switch for allocation passed into RTTOV subroutines

  ! Allocate input profile arrays
    ALLOCATE(profiles(nprof), stat = alloc_status(1))
    CALL rttov_error('mem allocation error for profile array' , lalloc = .TRUE.)

    CALL rttov_alloc_prof(     &
      errorstatus,             &
      nprof,                   &
      profiles,                &
      nlevels,                 &
      opts,                    &
      asw,                     &
      coefs = coefs(1),    &
      init = .TRUE._jplm)
    CALL rttov_error('allocation error for profile arrays')

! --------------------------------------------------------------------------
! 5. Store profile data in profile type
! --------------------------------------------------------------------------
    DO i = 1, NPROF
      profiles(i)%p(:) =  p_in(i, :)
      profiles(i)%t(:) =  t_in(i, :)
      profiles(i)%q(:) =  q_in(i, :)

      WHERE(profiles(i)%q(:) < 1e-4)
        profiles(i)%q(:) = 1e-4
      END WHERE

!      profiles(i)%o3(:) =  o3_in(i, :)
!      profiles(i)%co2(:) =  co2_in ! scalar
!      profiles(i)%ch4(:) =  ch4_in ! scalar
!      profiles(i)%n2o(:) =  n2o_in ! scalar
!      profiles(i)%co(:) =  co_in   ! scalar
           
      profiles(i)%cfraction  =  0._jprb
      profiles(i)%ctp        =  500._jprb
      
      ! 2m parameters 
      profiles(i)%s2m%p  =  p_surf(i)
      profiles(i)%s2m%t  =  t_in(i, 1) !t_surf(i)
      profiles(i)%s2m%q  =  q_in(i, 1) !q_surf(i)
      profiles(i)%s2m%u  =  2 !u_surf(i) ! DAR: hard-coded at 2ms-1?
      profiles(i)%s2m%v  =  2 !v_surf(i) ! DAR: hard-coded at 2ms-1?
      profiles(i)%s2m%wfetc  =  10000._jprb ! DAR: default?
      
      ! Skin variables for emissivity calculations
      profiles(i)%skin%t          =  t_skin(i)

      ! FASTEM coefficients - for MW calculations
      profiles(i)%skin%fastem(1)  =  3.0
      profiles(i)%skin%fastem(2)  =  5.0
      profiles(i)%skin%fastem(3)  =  15.0
      profiles(i)%skin%fastem(4)  =  0.1
      profiles(i)%skin%fastem(5)  =  0.3
    
      profiles(i)%zenangle      = zenang_in ! pass in from COSP

      profiles(i)%azangle       = 0._jprb ! hard-coded in RTTOV9 int

      profiles(i)%latitude      = latitude(i)
      profiles(i)%elevation     = h_surf(i)

      profiles(i)%sunzenangle   = 0._jprb ! hard-coded in RTTOV9 int
      profiles(i)%sunazangle    = 0._jprb ! hard-coded in RTTOV9 int

      ! surface type
      profiles(i)%skin%surftype  = lsmask(i)

      ! DAR: hard-coded to 1 (=ocean water) in RTTOV 9 int
      profiles(i)%skin%watertype = 1

      ! ! Cloud and aerosol profile data go here
      ! profiles(i)%aer_data   = .FALSE.
      ! profiles(i)%cld_data   = .FALSE.
      profiles(i) %idg         = 0._jprb
      profiles(i) %ish         = 0._jprb

      ! IF( lclouds ) THEN
      !   profiles(i) %cloud(:,:)  = 0._jprb
      !   profiles(i) %cfrac(:,:)  = 0._jprb
      ! ENDIF

      ! DAR: AEROSOL AND CLOUD PROFILE DATA GO HERE BUT NEED VALID COEFS    
    ENDDO
! End of 5.

    ich_temp = 1
    nchannels = nch_in
    DO no_id = 1, NRTTOVID

! --------------------------------------------------------------------------
! 3. Build the list of profile/channel indices in chanprof
! --------------------------------------------------------------------------

      ALLOCATE(nchan(nprof))     ! Number of channels per profile
      nchan(:) = SIZE(ichan(:,no_id))  ! = nch_in

    ! Size of chanprof array is total number of channels over all profiles
    ! Square in this case - here same channels done for all profiles
      nchanprof = SUM(nchan(:))
    
    ! Pack channels and input emissivity arrays
      ALLOCATE(chanprof(nchanprof))
!      ALLOCATE(emis(nchanprof))

      nch = 0_jpim
      DO j = 1, nprof
        DO jch = 1, nchan(j)
          nch = nch + 1_jpim
          chanprof(nch)%prof = j
          chanprof(nch)%chan = ichan(jch, no_id)
!          emis(nch)          = surfem(jch, no_id)
!         brdf(nch)          = input_brdf(jch)
          
!       ichan(:, 1)   =  ichan_in
!       surfem(:, 1)  =  surfem_in      
        ENDDO
      ENDDO
    
! End of 3.   

  ! Allocate output radiance arrays
      CALL rttov_alloc_rad( &
        errorstatus,        &
        nchanprof,          &
        radiance,           &
        nlevels - 1_jpim,   & ! nlayers
        asw)
      CALL rttov_error('allocation error for radiance arrays')
      
    ! Allocate transmittance structure
      CALL rttov_alloc_transmission( &
        errorstatus,             &
        transmission,            &
        nlevels - 1_jpim,        &
        nchanprof,               &
        asw,                     &
        init=.TRUE._jplm)
      CALL rttov_error('allocation error for transmission arrays')
      
    ! Allocate arrays for surface emissivity
      ALLOCATE(calcemis(nchanprof), stat=alloc_status(1))
      ALLOCATE(emissivity(nchanprof), stat=alloc_status(2))
      CALL rttov_error('mem allocation error for emissivity arrays' , lalloc = .TRUE.)
      
! DAR: where should this go?
!      surfem(:,:) =  0.0_JPRB
      emissivity(:)%emis_in = 0.0_JPRB       
!      input_emissivity(:) = emissivity(:)
      calcemis(:) = emissivity(:)%emis_in < 0.01_JPRB  
          
    ! DAR: This is a basic interface with RTTOV 11
    ! What's missing are: 
    !   optional arrays for calculating reflectances for the visible
    !   optional types for storing secondary radiances
    !                  and cloud and aerosol user-defined optical parameters
    !                  and PC calculations
    
    CALL rttov_direct(           &
      errorstatus,               &! out
      chanprof,                  &
      opts,                      &
      profiles,                  &! in
      coefs(no_id),              &! in
      transmission,              &! out
      radiance,                  &
      calcemis = calcemis,       &! in
      emissivity = emissivity)    ! inout

    CALL rttov_error('rttov_direct error')

    tbs(1:prof_num , ich_temp:ich_temp + SIZE(ichan(:,no_id)) - 1) = &
      TRANSPOSE(RESHAPE(radiance%bt(1:nchanprof), (/ SIZE(ichan(:,no_id)), prof_num/) ))

    ich_temp = ich_temp + SIZE(ichan(:,no_id))

  ! --------------------------------------------------------------------------
  ! 8. Deallocate all RTTOV arrays and structures
  ! --------------------------------------------------------------------------

!    DEALLOCATE (emis,        stat=alloc_status(1))
!   DEALLOCATE (brdf,        stat=alloc_status(2))
    DEALLOCATE (nchan,       stat=alloc_status(3))
    DEALLOCATE (chanprof,    stat=alloc_status(4))
    DEALLOCATE (emissivity,  stat=alloc_status(5))
    DEALLOCATE (calcemis,    stat=alloc_status(6))
!   DEALLOCATE (reflectance, stat=alloc_status(7))
!   DEALLOCATE (calcrefl,    stat=alloc_status(8))
    CALL rttov_error('RTTOV array deallocation error', lalloc = .TRUE.)

    asw = 0 ! Switch for deallocation passed into RTTOV subroutines

  ! Deallocate radiance arrays
    CALL rttov_alloc_rad(errorstatus, nchannels, radiance, nlevels - 1_jpim, asw)
    CALL rttov_error('radiance deallocation error')

  ! Deallocate transmission arrays
    CALL rttov_alloc_transmission(errorstatus, transmission, nlevels - 1_jpim, nchannels, asw)
    CALL rttov_error('transmission deallocation error')

  ENDDO

  ! Deallocate profile arrays
  CALL rttov_alloc_prof(errorstatus, nprof, profiles, nlevels, opts, asw)
  CALL rttov_error('profile deallocation error')

  DEALLOCATE(profiles, stat=alloc_status(1))
  CALL rttov_error('mem deallocation error for profile array', lalloc= .TRUE.)

  DO no_id = 1, nrttovid
    CALL rttov_dealloc_coefs(errorstatus, coefs(no_id))
    CALL rttov_error('coefs deallocation error')
  ENDDO

  DEALLOCATE(coefs, stat=alloc_status(1))
  CALL rttov_error('mem deallocation error for coef array', lalloc= .TRUE.)
    

  !======================================================================!
  ! 9.  Deallocate COSP-specific arrays                                  !
  !======================================================================!

  CONTAINS

    SUBROUTINE rttov_error(msg, lalloc)
      CHARACTER(*) :: msg
      LOGICAL, OPTIONAL :: lalloc

      IF(PRESENT(lalloc)) THEN
        IF(lalloc) THEN
          IF (ANY(alloc_status /= 0)) THEN
            WRITE(*,*) msg
            CALL rttov_exit(errorstatus)
          ENDIF
        ELSE
          IF (errorstatus /= errorstatus_success) THEN
            WRITE(*,*) msg
            CALL rttov_exit(errorstatus)
          ENDIF
        ENDIF
        
      ENDIF
    END SUBROUTINE rttov_error
  END SUBROUTINE RTTOV_MULTPROF
END MODULE MOD_COSP_RTTOV
